/*
9.27
1.function matrix_calc_taskA,matrix_calc_taskB
2."N" means iteration times
10.2
1.merge the for loop in taskB
3.move memory access variable out of "for loop".
10.4
1.add avx512 SIMD 
*/
#ifndef _CALC_H_
#define _CALC_H_

#ifdef __cplusplus
extern "C"{	
#endif

#include<omp.h>//multithread
#include<immintrin.h> //avx
#include<stdio.h>//printf
typedef struct {
	int* 			rowArray;
	const int* 		rowOffset;
	int 			rowArraySize;
	const int* 		columnIndice;
	const double* 	S;
	const double*	valueNormalMatrix;
	double* 		Id;
}TaskMatrixInfoA;

typedef struct {
	const double * 	valueSpiceMatrix;//N*2
	const int* 		rowOffset; 		//n+1
	const int* 		columnIndice; 	//N

	double* 		A;	//N
	double*		 	S;	//nx1
	double*			R;	//nx1
	double* 		H;	//nx1
	const double* 	D;	//N*2
	double* 		IC;	//nx1
	double* 		IG;	//nx1
	double 		alpha;
	int * 			rowArray; //kx1
	int 			rowArraySize;
	void* 			hdl;  
}TaskMatrixInfoB;

void matrix_calc_taskA(TaskMatrixInfoA** listDataList,int N)
{
//int coreNum = omp_get_num_procs();
//omp_set_num_threads(32);
//#pragma omp parallel
//	{
//#pragma omp for
	for (int i = 0 ;i < N; ++i)
	{	
		int row_array_size= listDataList[i]->rowArraySize;	
		TaskMatrixInfoA * list_Data_List=listDataList[i];
		const double *valueNormalMatrix = list_Data_List->valueNormalMatrix;
		const double *S = list_Data_List->S;
		const int *columnIndice = list_Data_List->columnIndice;

		
		for (unsigned int it = 0; it < row_array_size; ++it	)
		{
			const int node = list_Data_List->rowArray[it];
			const int column_num1 = list_Data_List->rowOffset[node];
			const int column_num2 = list_Data_List->rowOffset[node + 1];
			double &id=list_Data_List->Id[node];
			double id2 = id;
			//prepare for avx 
			__m512d m1,m2,m3,sum;			
			sum = _mm512_setzero_pd();
			 
			const int array_left = (column_num2 - column_num1)&0x07;
			int j;
			
			for (j = column_num1;j < column_num2 - array_left;j+=8)
			{
				m1 = _mm512_set_pd(
					valueNormalMatrix[j],
					valueNormalMatrix[j+1],
					valueNormalMatrix[j+2],
					valueNormalMatrix[j+3],
					valueNormalMatrix[j+4],
					valueNormalMatrix[j+5],
					valueNormalMatrix[j+6],
					valueNormalMatrix[j+7]
					);
				m2 = _mm512_set_pd(
					S[columnIndice[j]],
					S[columnIndice[j+1]],
					S[columnIndice[j+2]],
					S[columnIndice[j+3]],
					S[columnIndice[j+4]],
					S[columnIndice[j+5]],
					S[columnIndice[j+6]],
					S[columnIndice[j+7]]
					);
				m3 = _mm512_mul_pd(m1,m2);
				sum = _mm512_add_pd(sum,m3);
			}

			//多余的数据可以直接累加，不再用avx512寄存器 
			double tmp[8] __attribute__((aligned(64)));
			_mm512_store_pd(tmp,sum);
			id += tmp[0]+tmp[1]+tmp[2]+tmp[3]+tmp[4]+tmp[5]+tmp[6]+tmp[7];
			for(int ti = j;ti<column_num2;++ti){
				id +=valueNormalMatrix[ti]*S[columnIndice[ti]];
			}			
		}
	}
//	}
}

//Task B interface
void matrix_calc_taskB(TaskMatrixInfoB** listDataList,int N)
{
//int coreNum = omp_get_num_procs();
//omp_set_num_threads(32);
//#pragma omp parallel
//	{
//#pragma omp for
//	printf("\n\n************************************************\n\n\n");
	for (int i = 0;i<N;i++)
	{
		//指针变量
		TaskMatrixInfoB * list_Data_List = listDataList[i];
		int row_array_size = list_Data_List->rowArraySize;

		const double alpha = list_Data_List->alpha;
		double* A = list_Data_List->A;
		const double* S = list_Data_List->S;
		const int* columnIndice=list_Data_List->columnIndice;
		const double* valueSpiceMatrix=list_Data_List->valueSpiceMatrix;


		double 	tmp[8] __attribute__((aligned(64)));
		double 	ig;
		double 	ic;
		int 	p;
		int		u;
		int 	k;

		for(int it = 0;it < row_array_size;++it)
		{
			//指针变量
			int row = list_Data_List->rowArray[it];
			const int k1 = row << 1;
			double current = list_Data_List->D[k1];
			double charge = list_Data_List->D[k1+1];
			int column_num1 = list_Data_List->rowOffset[row];
			int column_num2 = list_Data_List->rowOffset[row + 1];
									
			double &IG = list_Data_List->IG[row];
			double &IC = list_Data_List->IC[row];
			int array_left= (column_num2 - column_num1)&0x07;

			p=column_num1;

	        __m512d COND,S_COL,CAP,ALPHA_VEC,TASK_4_5_7_8_IG,TASK_4_5_7_8_IC,TASK_9_A;
			ALPHA_VEC = _mm512_set1_pd(alpha);
			TASK_4_5_7_8_IG = _mm512_setzero_pd();
			TASK_4_5_7_8_IC = _mm512_setzero_pd();


			for (p = column_num1;p < column_num2-array_left;p+=8)
			{
				k = p << 1;
				// (4)(5) IG矩阵的乘法，即cond矩阵
				COND = _mm512_set_pd(
					valueSpiceMatrix[k+14],
					valueSpiceMatrix[k+12],
					valueSpiceMatrix[k+10],
					valueSpiceMatrix[k+8],
					valueSpiceMatrix[k+6],
					valueSpiceMatrix[k+4],
					valueSpiceMatrix[k+2],
					valueSpiceMatrix[k]
					);

				// (4)(5) S矩阵
				S_COL = _mm512_set_pd(
					S[columnIndice[p+7]],
					S[columnIndice[p+6]],
					S[columnIndice[p+5]],
					S[columnIndice[p+4]],
					S[columnIndice[p+3]],
					S[columnIndice[p+2]],
					S[columnIndice[p+1]],
					S[columnIndice[p]]
					);
				//IG[row]+=cond*S[col]
				TASK_4_5_7_8_IG = _mm512_add_pd(TASK_4_5_7_8_IG, _mm512_mul_pd(COND,S_COL));

				// (7)(8) IC矩阵的乘法，即cap矩阵
				CAP=_mm512_set_pd(
					valueSpiceMatrix[k+15],
					valueSpiceMatrix[k+13],
					valueSpiceMatrix[k+11],
					valueSpiceMatrix[k+9],
					valueSpiceMatrix[k+7],
					valueSpiceMatrix[k+5],
					valueSpiceMatrix[k+3],
					valueSpiceMatrix[k+1]
					);		
					
				// IC[row]+=cap*S[col]
				TASK_4_5_7_8_IC = _mm512_add_pd(TASK_4_5_7_8_IC, _mm512_mul_pd(CAP ,S_COL));
				//(9)不用512bit寄存器，直接算，减少来回加载
				A[p]  =valueSpiceMatrix[k]  +alpha*valueSpiceMatrix[k+1];
				A[p+1]=valueSpiceMatrix[k+2]+alpha*valueSpiceMatrix[k+3];
				A[p+2]=valueSpiceMatrix[k+4]+alpha*valueSpiceMatrix[k+5];
				A[p+3]=valueSpiceMatrix[k+6]+alpha*valueSpiceMatrix[k+7];
				A[p+4]=valueSpiceMatrix[k+8]+alpha*valueSpiceMatrix[k+9];
				A[p+5]=valueSpiceMatrix[k+10]+alpha*valueSpiceMatrix[k+11];
				A[p+6]=valueSpiceMatrix[k+12]+alpha*valueSpiceMatrix[k+13];
				A[p+7]=valueSpiceMatrix[k+14]+alpha*valueSpiceMatrix[k+15];
			}

			_mm512_store_pd(tmp,TASK_4_5_7_8_IG);
			IG += tmp[0]+tmp[1]+tmp[2]+tmp[3]+tmp[4]+tmp[5]+tmp[6]+tmp[7];
			current -= tmp[0]+tmp[1]+tmp[2]+tmp[3]+tmp[4]+tmp[5]+tmp[6]+tmp[7];

			_mm512_store_pd(tmp,TASK_4_5_7_8_IC);
			IC += tmp[0]+tmp[1]+tmp[2]+tmp[3]+tmp[4]+tmp[5]+tmp[6]+tmp[7];
			charge -= tmp[0]+tmp[1]+tmp[2]+tmp[3]+tmp[4]+tmp[5]+tmp[6]+tmp[7];
			//多余部分直接单个乘法，不用512bit寄存器 			
			for(u=p;u<column_num2;u++)
			{
				k = u << 1;
				ig = valueSpiceMatrix[k]*S[columnIndice[u]];
				ic = valueSpiceMatrix[k+1]*S[columnIndice[u]];
				A[u]=valueSpiceMatrix[k]+alpha*valueSpiceMatrix[k+1];
				IG +=ig;
				IC +=ic;
				current -= ig;
				charge -= ic;
			}
			list_Data_List->R[row] = current; 
			list_Data_List->H[row] = charge;
		} 
	}
//	}
}

#ifdef __cplusplus
}
#endif

#endif
